function y = computeDeltaLambda(element,normxi,alpha)
dlambda = 0;
a = alpha;
tol = 1e-6;
mu = 0.5*element.E/(1+element.nu);
while 1
    H1 = computeKinematicHardeningModulus(element,alpha);
    H2 = computeKinematicHardeningModulus(element,a);
    g = -sqrt(2/3)*computeIsotropicHardeningModulus(element,a)...
        + normxi - (2*mu*dlambda+sqrt(2/3)*(H2-H1));
    if abs(g) < tol
        break;
    end
    dH = computeKinematicHardeningSlope(element,a);
    dK = computeIsotropicHardeningSlope(element,a);
    dg = -2*mu*(1+(dH+dK)/(3*mu));
    dlambda = dlambda - g/dg;
    a = alpha + sqrt(2/3)*dlambda;
end
y = dlambda;
end
